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Abstract 

We present an analysis of the stability, energy and torque properties of 
a model Bursian diode in a one dimensional Eulerian framework using the 
cold Euler-Poisson fluid equations. In regions of parameter space where there 
are two sets of equilibrium solutions for the same boundary conditions, one 
solution is found to be stable and the other unstable to linear perturbations. 
Following the linearly unstable solutions into the non-linear regime, we find 
they relax to the stable equilibrium. A description of this process in terms 
of kinetic, potential and boundary-flux energies is given, and the relation to 
a Hamiltonian formulation is commented upon. A non-local torque integral 
theorem, relating the prescribed boundary data to the average current in the 
domain, is also provided. These results should prove useful for understanding 
Bursian diodes in general, as well as for control applications and benchmarking 
numerical codes. 

1 Introduction 

In its simplest form, a diode consists of two conducting electrodes with a relative elec- 
tric potential bias \<f>i\, and a distribution of moving charge carriers. Fundamentally, 
the transport of these charge carriers is constrained, self-consistently, by non-linear 
space charge effects. For example, in the case of a steady un-neutralized electron flow 
in one dimension (a Bursian diode), the charge flux cannot exceed the analytically 
derivable 'Child-Langmuir limit' that depends only on \<fii\, the size of the domain 
and the velocity of the incoming electrons [H El E] . Mechanistically, if the electron 
flux exceeds the limiting value, there is a charge build-up - a virtual cathode - and 
an associated electric field that resists the passage of additional electrons. 

Understanding and controlling the onset of this virtual cathode, as well as other, 
nearby, physically and numerically accessible states, their stability properties, and 
the energy demands of maintaining a diode flow, has applications in a wide range of 
settings that are well reviewed by Ender et al. [I]. Some examples include inertial- 
electrostatic confinement [5]; pinch reflex diodes for intense ion beam generation [6]; 
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vircators [7J; reflex triodes for microwave generation [8]; photoinjectors [SI [10], and 
producing GHz to THz electromagnetic radiation [TT| [12]. 

Historically, much of the illuminating analysis has come from simulations, espe- 
cially in complex geometries and for kinetic systems. To ensure the fidelity of future 
codes, a good understanding of the basic physics and a suite of test cases for bench- 
marking is desirable. Furthermore, in the fluid limit, diodes are readily analyzable, 
energetically open system, as the entering and exiting particles carry with them ki- 
netic and potential energy. They therefore constitute a good practical example from 
which to extend the Hamiltonian description of a fluid beyond energetically closed 
systems [T3[ [II] . To these ends, this paper investigates the linear and non-linear 
dynamics, and the time dependent energy evolution of the two-equilibria region of 
parameter space supported by the simple Bursian diode above. 

The plan of this paper is as follows. In section [2] we introduce our equations and 
review what is known about their time-independent, i.e. equilibrium, solutions. We 
focus on regions of parameter space that supports two distinct equilibria. In section 
[3] we present a new perspective on their linear stability, showing one to be stable 
and the other unstable. In section [4] we continue our investigation by following the 
linear instability into the non-linear regime, and discuss the associated system energy 
and torque, and their role as diagnostics. In section [5] we conclude and discuss some 
potential applications for our results. 

2 Equilibrium solutions 

A standard model to describe charge carriers in an electric potential is the cold ID 
hyperbolic-elliptic Euler-Poisson system given by 

d tP + d x {pv) = 0, (1) 
d t v + vd x v = d x <p, (2) 
d xx <P = P, (3) 

where x,t,p,v,<p are the scaled position, time, density, velocity and potential for an 
electron fluid and the time independent Dirichlet boundary conditions are 

p(x = 0) = po, v(x = 0) = v , cf)(x = 0) = 0, <p(x = d) = (pi. (4) 

We normalize using 

(x,t,v,<f>,p) = (x'/L,t'/T,v'/(L/T)^'/ V ,p'/R) 
<p = (m e /q e )(L/T) 2 R = (e /q e )($/ L 2 ). 
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Figure 1: Potential profiles associated with branch I (solid lines) and branch II 
(dashed lines) equilibria for (d,i> , 0Oj 0i) = (1.0,1.0,0,0.25), and po = 2.45 (red, no 
symbols), po = 2.00 (blue, diamonds), and po = 1-00 (black, squares) respectively. 



The primed variables are unsealed, L,T are characteristic length and time scales, m e 
is the electron mass, q e the fundamental charge (positive), e is vacuum permittivity, 
and the electric field is —d x cj). 

In the steady state ([T| and ([2| constrain the current pv and the energy density 
of a fluid element, kinetic plus potential v 2 /2 — <p, to be constant across the domain 
(the minus is because electrons are negatively charged). This implies that for given 
boundary conditions i.e. (HI), the two unspecified fields at the outgoing boundary 
p(d),v(d) are uniquely determined. The motion of the fluid can be understood ener- 
getically in terms of Hamilton's principle, the principle of least action, from which ^ 
can be derived. The gain (loss) in the kinetic energy of a fluid element as it crosses 
the domain equals its loss (gain) in potential energy as work is done on (against) 
it by the electric field (that accelerates electrons from the emitting cathode to the 
collecting anode, in the case of a monotonically increasing potential). 

It is known that there are two kinds of equilibrium solutions to these equations 
[3], and we review them here. For <pi > 0, their implicit expressions are given by: 
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(5) 
(6) 



where $ = \/l + 2<P/vq and $d = y/l + 2<Pi/vq are normalized potentials, and 
6 = 4/3(1 + $;f ) < 6 = 4/3($ d + 2)($ d - l) 1 ' 2 <& = 4/3 (1 + $ d ) 3/2 demarcate 
(non-exclusively) the boundaries between solutions monotonic in <p given by ^ and 
solutions with a single turning point given by ([6]). 
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To close these equations, a is needed. It satisfies 



($ d - 2a) v^d + a = - A ^^d ± (1 - 2a) + a, (7) 

where the positive and negative signs correspond to ^ and (|6]) respectively. 

Necessary and sufficient conditions for determining the number of solutions, zero, 
one or two, are given succinctly by 

d\J 8p / w o > £3 : zero solutions, (8) 

d\J Sp Q /vQ < £1 or dyj 8p /vQ = £ 3 : one solution, (9) 

£1 < 8po/vQ < £3 : two solutions. (10) 

The number of accessible solutions is a function of d, v , p , and (pi, Figs. [I] and 
[2j For example, for (d, t>o,0i) = (1,1,0.25), there are no steady state solution for 
Po > 2.5, two when 2.5 > po > 1-2, and one when p < 1-2. We denote the solution 
with larger as branch I, the other as branch II. In the literature, these are known 
as the C-branch and C-overlap branch respectively [T5] . 

It is the stability, dynamics and energy of perturbations to the equilibria in the 



region of parameter space given by (10), that are the focus of this paper. While 
these have been investigated before in a Lagrangian framework, our approach in an 
Eulerian framework is new, and has several advantages. Specifically, it allows for a 
direct interpretation of solutions that are functions of x and t, rather than Lagrangian 
coordinates; the discrete nature of the linear eigenmodes are a natural product of the 
formulation; and the description is robust to changes that would not allow for a 
Lagrangian analysis. 

In accordance with earlier studies, we find that the C-overlap branch is unstable 
to linear perturbations, and we follow these into non-linear regime [3j [T71 [18] . 



3 Linear stability analysis 

We wish to determine the linear stability properties of branch I and branch II equi- 
librium solutions to @-((3]) when (10) holds. To do so, we use ([5])-([7| subject to 



^j) to construct equilibria p(x), v(x), (f)(x), and to these we add small perturbations 
(8p,5v,5(f>) = (Sp(x),Sv(x),S4>(x))e xt that obey (8p,8v,8(f)) = (0,0,0) at x = and 
86 — at x — d. 
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Figure 2: The number of equilibrium solutions to (JTh-QSJ) : Zero - horizontal lines, one 
- vertical lines, two - crossing lines, for varying Po,v ,(f)i. Left: (d,po) = (1, 1) Right: 
(d, vq) = (1, 1). The zero and two solution boundary - d^JSpo/vl — ^ 3 - corresponds 
to the space charge limiting current derived by Child and Langmuir for vq = 0, Jaffe 
for vq > 0, and recently, for time varying solutions, by Caflisch and Rosin [T| |2| 131 fl6] . 
The circled dot corresponds to the parameters used in Fig. [3] and the star to those 
used in Figs. [5] and [6] 



Linearizing, we obtain 



X5p + d x (p5v) + d x (v5p) = 0, 

X5v + d x (v8v) = d x 5(p, 
d xx S(j) = 5p. 



(11) 
(12) 

(13) 



This system can be written as 



where 



A :-- 



8p 
6v 



A ■ 



-d x v — vd x 
d x {d xx )- x 



8p 
Sv 



-d x p - pd x 
-d x v — vd x 



(14) 



(15) 



The eigenvalues of (14), determine the linear stability of the system, 3ft(A) > de- 



scribes unstable modes, and 5f (A) < 0, stable modes. To compute A, we discretize the 



operator matrix (15) using three methods: a uniform grid with an upwind scheme; a 
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Figure 3: Left: Eigenvalues associated with perturbations to branch II equilibria 
with parameters (d, po, Vo, <pi) = (1,2.4,1,0.25). Only a single positive eigenvalue, 
the first one, exists, corresponding to an unstable, purely growing mode. The re- 
maining eigenvalues are in complex conjugate pairs with 9ft(A) < 0, corresponding to 
damped, traveling waves. Results are calculated using Chebyshev spectral methods 
with N = 100, 200, 400 modes corresponding to black circles, triangles, and diamonds 
respectively - which overlap completely. 

uniform grid with a centered difference scheme; and a Chebyshev grid with an asso- 
ciated polynomial interpolation [TH] . The discrete spectrum of eigenvalue-eigenvector 
solutions - a discreteness not generally emphasized in the dispersion relations arising 
from Lagrangian analyses e.g. [TT1 ITS] . - are obtained numerically and shown in 
Figures |3j [4] and |5j All three schemes converge to the same result. 

Conducting a parameter scan, for branch II we find A > G 3ft for the first 
eigenvalue, the one with a single zero in the corresponding eigenfuctions. For the 
remaining eigenvalues in branch II, and all of branch I, 9ft(A) < 0. The system 
supports a single unstable mode. For example, for (d, po, Vq, <f>i) = (1, 1.5, 1,0.2), the 
most positive eigenvalues from branch II and I are 1.1 and —2.1 respectively - one 
mode is unstable, and the other stable. As the two equilibrium solutions merge at 
d\/8po/ v o — ^3; the unstable eigenvalue of branch II obeys 9ft(A) — > 0. Approaching 
the other boundary of the two solution region d^/8p /vQ = £i, the full-width, half- 
maximum of the corresponding eigenmode — > 0. It remains to be seen whether this 
singular mode bears any fundamental relation to the singularity that forms in the 
case that the current exceeds the Child-Langmuir limit [201 [16] . 

4 Nonlinear dynamics 

For small time, coupling between infinitesimal amplitude perturbations, and their 
feedback on the equilibrium solutions, is negligible. However, because A > for 
one mode, that mode grows exponentially and the perturbations quickly reach non- 
linear amplitudes. In this case, the methods and results of section [3] are no longer 
applicable. In the non-linear regime, the most general method for solving ([l|-((3]) 
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Figure 4: The most positive eigenvalues of (14) i.e. A, associated with perturbations 
to branch II (9?(A) > 0, unstable) and branch I (3?(A) < 0, stable) solutions for the 
parameters (vq, d) = (1, 1), and varied po and <pi - see Fig. [2j 



is numerical integration; although the method of characteristics can also be used to 
obtain complete analytic solutions in a Lagrangian framework [TBI II]- The method 
used here, an Eulerian approach, has the advantage that it is naturally formulated 
as a two point Dirichlet boundary value problem for 0, which can easily be realized 
experimentally. The alternative Lagrangian approach is more naturally formulated 
as a Cauchy problem including d x <p, which is harder to realize experimentally, and 
from which the corresponding Dirichlet conditions are non-trivial to obtain [16J. 

We favor the numerical approach. We employ MacCormack's method to integrate 
the hyperbolic equations Q-Q, and solve the elliptic Poisson equation ^ at each 
time step using a finite-difference description and inverting a tridiagonal matrix. Our 
simulations are initialized with unstable equilibrium solutions from branch II and 
numerical noise provides broadband perturbations which are constrained to obey Q . 
The solutions to our perturbed system are well matched by our linear results for 
small time, and in the final state the solutions have relaxed to the stable branch 
I equilibrium solutions with the same boundary conditions as the initial, unstable 
equilibrium, Fig. [5j 

Physical insight and a set of numerical benchmarks for the system, can be obtained 
by considering both the energetics of the system and its global torque. In the next 
section, we examine each in turn, and derive a set of integral equations that describe 
the system's spatially averaged properties and their interaction with the boundaries. 

These type of equations are both less computationally demanding to solve (which 
is unimportant here, but may matter in higher dimensions or kinetic models), and 
do not require knowledge of the fundamental unaveraged solutions. Furthermore, 
being able to relate prescribed boundary value data to derived domain data offers a 
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Figure 5: Left hand side: Evolution of perturbed field quantities associated with 
the unstable branch II equilibrium for t < 5.6 and steady boundary conditions 
(d, po, v , 0i) = (1, 3, 1, 1). Snapshots are every t = 0.35 starting from the 8p, 5v, 5<p = 
initial conditions. The unstable linear eigenmodes (red squares) with growth rate 
exp(At) - section [3]- match the fully non-linear solutions for small time. Right hand 
side: Evolution of full solutions to ([]])- ([3]) starting at the same branch II equilibrium, 
for t < 11.5. Snapshots are every t = 0.5, and the initial state (red circles) is given 
by © - ([7]). The final state (blue diamonds) is the same stable branch I equilibrium 
derived from the same set of equations and boundary conditions as the initial state. 



new avenue for both control, and experimental measurement of spatially-distributed 
system properties. 



4.1 Energetics 

Even at the model equation level considered here, energy insights may be important 
for industrial purposes j7]. In this section, we examine the evolution and balance of 
the standard energy integrals. We leave further detailed discussion to a forthcoming 
paper in which we present a tailored Bursian diode-battery model |21j . 

We start by multiplying ^ by v and combining it with to yield an evolution 
equation for the kinetic energy JC = pv 2 /2 balance of the system 

dtK + d x (vJC) = pvd x( j), (16) 

where pvd x (p is the negative Joule heating term. Integrating over x, the total kinetic 
energy in the system is given by 

d t K, = v K, - v d IC d + pvd x (j) (17) 
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Figure 6: Evolution of the total (spatially integrated) energy of a Bursian diode 
system from an initial unstable equilibrium solution (branch I) to a final stable equi- 
librium solution (branch II) for the same set of boundary conditions as in Fig. [5] 
Whilst the final total energy state £ is slightly greater than the initial state, the 
Hamiltonian % is a strictly decreasing function of time. The system's dominant form 
of energy switches to kinetic K, from internal potential energy Vj as time progresses. 



where over-bars denotes spatially integrated quantities J Q d dx and subscripts 0, d indi- 
cate that the associated quantity is to be evaluated at x — 0, d respectively. There are 
two contributions to the total kinetic energy: the difference in the boundary fluxes of 
kinetic energy, and the work done on the fluid by the electric field. 

To describe the total energy balance in the system, it helps to decompose = 0e + 
4>\ into external and internal components, and these satisfy Laplace's and Poisson's 
equations respectively: 



d xx <pE = 0, with E (O) = 0, 4> E (d) 
d xx <pi = p, with 4>i (0) = 0, 4>i(d) = 



0. 



The solution to (18) is simply (f>E 



i/d)x, and the Green's function for <pi is 



dx'p(x') 



o 



XX 



Rewriting (16) in conservative form using ([T]) and (19), we have 

d t S + d x (v (8 + Vi)) = 0, 



(18) 
(19) 



(20) 



(21) 



where £ = K + Ve + V\ is the combined kinetic /C, external potential Ve = —p4>E 
and internal potential V\ = — p0i/2 energy of a fluid element, and we have made use 
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of the fact that 0e,0i(O) and 4>i(d) are time independent. Physically, the factor of a 
half in the definition of V\ is to avoid double counting particle interaction energies 
[22] • Mathematically, it arises from the symmetry properties of the Green's function 



(20) under x <^ x'. 



In the absence of net boundary fluxes, the second term in (21) vanishes upon 
integration. In this case, the total energy £ is conserved, and coincides with the fluid 
Hamiltonian H = pv 2 /2 + (d x (p) 2 /2, from which the equation of motion (2| can be 
derived [2"3T ITT] . The evolution of the various energy quantities is plotted in Fig. [6] 

Considerable work has been done on the non-linear stability of closed plasma and 
fluid systems using variational principles e.g. [2H ?, [25l [26]. However, for open 
systems i.e. ones with sources, like boundary fluxes, stability proofs are difficult to 
construct, and we do not attempt to do so here. Nevertheless, the non-linear stability 
and Hamiltonian structure of such systems has been the focus of recent work, and so 
a theorem tailored to the problem described here may be forthcoming [23 12HI 123 ED] • 

4.2 Torque and boundary conditions. 

Unlike energy, torque is not generally considered as an important property of diode 
systems. However, it is frequently invoked in describing stellar systems governed by 
([I]) - (J3|, in the context of which ^ is known as the Jeans equation, and <\> is the 
gravitational potential. We consider it here too and derive a simplified lower moment 
analogue to the astrophysical virial theorem including boundary effects [31J. As for 
the virial theorem, we find a 'basic structural relation that the system must obey' 



To proceed, we note that, uniquely, the ID version of Poisson's equation ^ can 
be directly integrated to yield 

PX 

dxd xx <f> = d x( f>{x) - d x <f>(0) = / dx'p(x') = M(x):=M x , (22) 



which is the mass between and x (which can vary with time). It follows from (22) 
that 

<h. = J dx (^J g dx'p(x') + d x <p(0)^J =M d {d-x) + d x <p(Q)d, (23) 
where x = dxxpj dxp is, by definition, the center of mass, and M d is defined in 

ph. 



Equation (23) has a very simple interpretation. By definition, the torque about a 
point d is T = Fr where r is the magnitude of the directional vector joining d and the 
point at which F, the force perpendicular to this vector, acts. We consider a force 
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acting at and proportional to the system's center of mass M^, and perpendicular to 
v x , say a gravitational force F = M d g. In this case, we have T = M d g(d — x), and 
so 



<h-d x <l>(a)d= (d-x)M d = T, 



(24) 
■d x <KQ)d, 



where we have absorbed g into the definition of T. For time independent 0i • 
this implies that the total torque on the system is constant. 

This results in an interesting relation between the current, the rate of change of 
the incoming electric field dtd x <f)(0) and exiting potential d%<$>\. Differentiating (24), 



d t T = -d t xM d + (d-x) d t M d = d t [fa - d x <t>{Q)d\ 
and, using Q, we have 



d t M d = - dxd x (pu) = pqUq - p d u d , 
Jo 

d t xM d = -dp d u d + J(d) - d t M d x, 



(25) 

(26) 
(27) 



where J(d) = J Q dx pu is the current in the domain, and (26) simply states that the 
rate of change of mass is the flux in minus the flux out. 



Combining (25), (26) and (27) we find 



PqUq = d (J(d) + d t <pi) - d t d x (j)(0), 



(2f 



which is the main result of this section^ 

Equation ( [28) ) relates the average current in the domain, a derived quantity, to a set 
of boundary data. This, potentially, affords a new avenue for control. As mentioned 
earlier, because ([T]) - ^ can be written in characteristic form, in a mathematical 
sense, the appearance of the incoming electric field — d x <p(0) is a more natural choice 
of boundary condition than 0i. 



5 Conclusion 

While Bursian diodes have been well studied over the last century, the advent of large 
scale, multi-dimensional particle in cell codes, and fluid codes in complex geometries 
have the potential to offer new insights into their basic physics and to guide their 

^^An alternative derivation of this result due to R.Caflisch can be obtained by simply evaluating 
the Green's function solution of ^ for cf> at x = d with appropriate boundary conditions - private 
correspondence. 
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design. The work provided here, whilst relatively basic in that it is one dimensional 
and uses a minimal set of equations, is thorough and, as such, provides a reliable 
set of results and integral theorems against which the results of simulations of more 
complicated systems can be compared. Where we have employed numerical tools we 
have cross checked our results using multiple methods and conducted appropriate 
convergence studies. 

Our results include a linear stability analysis of the unstable branch II equilibrium 
(C-overlap branch), and non- linear simulations of its evolution. We have found that 
its relaxed state is that of the stable branch I equilibrium with the same boundary 
conditions. We have also provided a quantitative discussion of the role of energy and 
torque in diagnosing and controlling the system, and, to the best of our knowledge, 
our interpretation of the latter is new in the literature. 

Possible extensions to this work include constructing a non-linear stability theorem 
in the spirit of Bernstein et al., but including boundary fluxes [21]; using the results 
herein for benchmarking more complicated systems including investigating the sta- 
bility of Child-Langmuir limited solutions; and prescribing optimizing and efficiency 
enhancing conditions or frameworks for diode operation [331 [16] . 
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